{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 5000000/5000000 [00:59<00:00, 83837.01it/s]\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "from tqdm import tqdm\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "NUMBER_OF_STEPS = 5000000\n",
    "SAMPLE_INTERVAL = 10\n",
    "\n",
    "length = int(NUMBER_OF_STEPS/SAMPLE_INTERVAL)\n",
    "x_list = np.zeros(length, dtype=float)\n",
    "v_list = np.zeros(length, dtype=float)\n",
    "inv_list = np.zeros(length, dtype=float)\n",
    "\n",
    "def nhc(M, pos_eta, vel_eta, mas_eta, Ek2, kT, dN, dt2):\n",
    "    dt4 = dt2 * 0.5\n",
    "    dt8 = dt4 * 0.5\n",
    "\n",
    "    # Update velocity of the last (M - 1) thermostat:\n",
    "    G = vel_eta[M - 2] * vel_eta[M - 2] / mas_eta[M - 2] - kT\n",
    "    vel_eta[M - 1] += dt4 * G\n",
    "\n",
    "    # Update thermostat velocities from M - 2 to 0:\n",
    "    for m in range(M - 2, -1, -1):\n",
    "        tmp = math.exp(-dt8 * vel_eta[m + 1] / mas_eta[m + 1])\n",
    "        G = vel_eta[m - 1] * vel_eta[m - 1] / mas_eta[m - 1] - kT\n",
    "        if m == 0:\n",
    "            G = Ek2 - dN * kT\n",
    "        vel_eta[m] = tmp * (tmp * vel_eta[m] + dt4 * G)\n",
    "\n",
    "    # Update thermostat positions from M - 1 to 0:\n",
    "    for m in range(M - 1, -1, -1):\n",
    "        pos_eta[m] += dt2 * vel_eta[m] / mas_eta[m]\n",
    "\n",
    "    # Compute the scale factor\n",
    "    factor = math.exp(-dt2 * vel_eta[0] / mas_eta[0])\n",
    "\n",
    "    # Update thermostat velocities from 0 to M - 2:\n",
    "    for m in range(M - 1):\n",
    "        tmp = math.exp(-dt8 * vel_eta[m + 1] / mas_eta[m + 1])\n",
    "        G = vel_eta[m - 1] * vel_eta[m - 1] / mas_eta[m - 1] - kT\n",
    "        if m == 0:\n",
    "            G = Ek2 * factor * factor - dN * kT\n",
    "        vel_eta[m] = tmp * (tmp * vel_eta[m] + dt4 * G)\n",
    "\n",
    "    # Update velocity of the last (M - 1) thermostat:\n",
    "    G = vel_eta[M - 2] * vel_eta[M - 2] / mas_eta[M - 2] - kT\n",
    "    vel_eta[M - 1] += dt4 * G\n",
    "\n",
    "    return factor\n",
    "\n",
    "\n",
    "# Main simulation loop\n",
    "dN = 1.0\n",
    "dt = 0.01\n",
    "dt2 = dt * 0.5\n",
    "kT = 1.0\n",
    "x = 0\n",
    "v = 1\n",
    "mass = 1.0\n",
    "k_spring = 1.0\n",
    "f = -k_spring * x\n",
    "M = 4\n",
    "tau = dt * 100\n",
    "pos_eta = [0.0] * M\n",
    "vel_eta = [0.0] * M\n",
    "mas_eta = [0.0] * M\n",
    "\n",
    "vel_eta[0] = vel_eta[2] = +1.0\n",
    "vel_eta[1] = vel_eta[3] = -1.0\n",
    "\n",
    "for i in range(M):\n",
    "    pos_eta[i] = 0.0\n",
    "    mas_eta[i] = kT * tau * tau\n",
    "    if i == 0:\n",
    "        mas_eta[i] = dN * kT * tau * tau\n",
    "\n",
    "Ek2 = 0.0\n",
    "factor = 0.0\n",
    "\n",
    "\n",
    "for step in tqdm(range(NUMBER_OF_STEPS)):\n",
    "    inv = mass * v * v * 0.5 + k_spring * x * x * 0.5\n",
    "    inv += kT * dN * pos_eta[0]\n",
    "\n",
    "    for m in range(1, M):\n",
    "        inv += kT * pos_eta[m]\n",
    "\n",
    "    for m in range(M):\n",
    "        inv += 0.5 * vel_eta[m] * vel_eta[m] / mas_eta[m]\n",
    "\n",
    "    if step % SAMPLE_INTERVAL == 0:\n",
    "        nn = int(step/SAMPLE_INTERVAL)\n",
    "        x_list[nn] = x\n",
    "        v_list[nn] = v\n",
    "        inv_list[nn] = inv\n",
    "\n",
    "    Ek2 = v * v * mass\n",
    "    factor = nhc(M, pos_eta, vel_eta, mas_eta, Ek2, kT, dN, dt2)\n",
    "    v *= factor\n",
    "\n",
    "    v += dt2 * (f / mass)\n",
    "    x += dt * v\n",
    "    f = -k_spring * x\n",
    "    v += dt2 * (f / mass)\n",
    "\n",
    "    Ek2 = v * v * mass\n",
    "    factor = nhc(M, pos_eta, vel_eta, mas_eta, Ek2, kT, dN, dt2)\n",
    "    v *= factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAokAAADUCAYAAADwWLkhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7SUlEQVR4nO3deXhU5dnH8e+dyZ6QsAUQQgj7vggBFBRkE1AQd3GvuKFWq9aq2Ne6ttpqrS1SFXGhLkUrqCBYJIAgKEoA2VcBIbJGdhKyPu8fCTRkmQnJzDxnZu7Pdc1FZs6Zc35JuJMnZ7kfMcaglFJKKaVUaWG2AyillFJKKefRQaJSSimllCpHB4lKKaWUUqocHSQqpZRSSqlydJColFJKKaXK0UGiUkoppZQqx/ogUURcIrJCRD63nUWpQCEiw0Rko4hsEZFHK1h+gYgcFpEfSh5/sJFTqUClNaYUhNsOAPwGWA8k2A6iVCAQERcwARgCZAJLRWS6MWZdmVW/NsaM8HtApQKc1phSxaweSRSRZOBiYJLNHEoFmF7AFmPMVmNMHjAFGGU5k1LBRGtMKeyfbn4ZeBgospxDqUDSBNhZ6nlmyWtlnSsiK0XkCxHp6J9oSgUFrTGlsHi6WURGAPuMMctE5AI3690B3AEQFxfXo23bdv4JqJQby5cvyzLGJFnavVTwWtn5NZcDzYwxx0TkIuBToHWFGytTY+3aaY0p+5YtC44a0/pSTlTV+rJ5TWJf4JKS4ooGEkTkPWPMDaVXMsZMBCYC9OiRZhZ/l+H/pEqVERMhP1ncfSbQtNTzZGBX6RWMMUdKfTxLRP4pIvWNMVllN1a6xtLS0kxGhtaYsk8kOGpM60s5UVXry9rpZmPMOGNMsjEmFRgNzCs7QFRKVWgp0FpEmotIJMX1M730CiLSSESk5ONeFNf6L35PqlRg0hpTCmfc3ayUOgPGmAIR+TUwG3ABbxlj1orI2JLlrwFXAneJSAGQA4w2xpQ9XaaUqoDWmFLFHDFINMZ8BXxlOYZSAcMYMwuYVea110p9/Arwir9zKRUstMaUsn93s1JKKaWUciAdJCqllFJKqXJ0kKiUUkoppcrRQaJSSimllCpHB4lKKaWUUqocHSQqpZRSSqlydJColFJKKaXK0UGiUkoppZQqRweJSimllFKqHB0kKqWUUkqpchwxLZ9yz9NsoAZ704UK4n65+8VKKaWUcig9kqiUUkoppcrRQaJSSimllCpHB4lKKaWUUqocHSQqpZRSSqlydJColFJKKaXK0UGiUkoppZQqRweJSimllFKqHO2T6CdF7podemhzWFTTPok1aKMoHhodirjfuLt3a49FpZRSyrn0SKJSSimllCpHB4lKKaWUUqocHSQqpZRSSqlydJColFJKKaXK0UGiUkoppZQqRweJSgUgERkmIhtFZIuIPOpmvZ4iUigiV/ozn1KBTmtMKR0kKhVwRMQFTACGAx2Aa0WkQyXr/RmY7d+ESgU2rTGlimmfxCpy1+YQPPcqLCqqfFmhh0aInpYXuNs47rN76kXoCnO/Qk2Wuzz2QdQ+ipXoBWwxxmwFEJEpwChgXZn17gWmAj39Gy90FRQUkJeXR2xsrO0oqma0xhwqJyeHiIgIwsN1+OIPeiRRqcDTBNhZ6nlmyWuniEgT4DLgNU8bE5E7RCRDRDL279/v1aCh5A9/+APR0dHExcXRvXt3pk2bZjuSqj6v1ZjWl3fk5ubSqVMnYmNjSUhIYPTo0WzZssV2rKCng0SlAk9Fx1DLHi9+GXjEGFPoaWPGmInGmDRjTFpSUpI38oWkPn36cP/99/Pkk0+Sm5vLFVdcwf3334/xdBpCOZHXakzryzuioqIYNWoUzz77LLfccgszZ86kR48ezJ0713a0oKbHa5UKPJlA01LPk4FdZdZJA6aUTKtYH7hIRAqMMZ/6JWGImDlzJvPmzePFF19k2LBhDBs2DID/+7//4+GHH6Z9+/Yep7ZUjqQ15gDGGO68806uv/56+vfvzx//+MdTyx555BHuvfdemjZt6mYLqqZ0kKhU4FkKtBaR5sDPwGjgutIrGGOan/xYRN4BPtdfXt61Y8cObrzxRlJTU8nNzeWqNzJOXyH1MjblwGfjF1FUkM/MBwbYCaqqQ2vMAV5//XXeeOMNmjdvTv/+/U9blpKSwmeffXbqeV5eHpGRkf6OGPT0dLNSAcYYUwD8muI7KtcDHxlj1orIWBEZazddaDDGcOutt5Kfn8+HH35IdHR0pevuWb2Ir/54PQcPHvRjQlUTWmP2bd26lQceeIChQ4fyyCOPVLqeMYYxY8Ywdqx+W3zB2pFEEWkK/AtoBBQBE40xf7eVR6lAYoyZBcwq81qFF9AbY37lj0yh5NNPPyU9PZ3x48fTunVrt+vG1GlI9oG9PPXUU7z88sv+CahqTGvMjpHjFwGQMen35BcJBX3vZNSEb05bZ8a95536WESoX78+L7zwAnfddRc9e+qN5t5k83RzAfBbY8xyEakFLBOROcaYsi0GlFLKMYqKivjd735Hp06dqnT0IjG5NSl9RvLKK69w11130bZtWz+kVCpwHdi6ij2rFtJuxB3E1C5/s8/JgeRJ+fUHEVXrTYZeM4Y+9//z1HXApQeTqnqsDRKNMbuB3SUfHxWR9RS3GLAySKxpH0RPvQzzCypfnlfovs/hiTz3N6jmFbh/v7toHtocEhnu/oqE6AiX2+VREW7e7/6tnvsoGu2jqPwvLCyMGTNmkJOTU+VebW0vupWsFXP405/+xOTJk32cUKnAVie1Iz1u/SMNOvSu0voRMXG0GX4Lqz/6K1mblpHUNs3HCUOHI65JFJFU4GzguwqW/a/HVJb2mFJK2de+fXu6d+9e5fWjatXhzjvv5P3332fXrrI3ySqlSpMwF2d17YcrIqrK70nufRFRifXZOvffPkwWeqzf3Swi8RR3rL/fGHOk7HJjzERgIkCPHmnacEwp5XcnT2/tWbWQzKVf0mX074iMSzyjbTz88MNcddVVNG7c2BcRlQoKKz94nlpnNafFgGvO6H2uiEh63PIUcUnaEsebrB5JFJEIigeI7xtjdHoCpZSjbZ3/EUcyNxMRE3/G723UqBHnnnuuD1IpFRy2bdvGzu9mkZ99tFrvr9uiC1G16ng5VWizNkiU4itL3wTWG2NespVDKaWq4uie7Rz4cSXNzhuFhHm4oLYSOTk53HbbbXpdolIVmDRpEiCk9Lmk2ts4+NM6lrxyP3nHy52YVNVg80hiX+BGYKCI/FDyuMhiHqWUqtSOb2YgrnCSew2v9jaio6PJyMjgb3/7m07Xp1Qp+fn5vPXWWzTocA4xdRpUezuu8EiyNi0j87tZnldWHlkbJBpjFhljxBjTxRjTreSh31WllOMU5ueSuXQ2jTqfX6PTWSLC2LFjWblyJUuXLvViQqUC28yZM9mzZw/NanAUESChSSvqNO/ET9/M0D/EvMD6jStO4csWNwAn8itvY3MoO9/tew97WH4sv8DtcnctdiI89JmpHeV+mqPE2Ai3y+OKKv8vFhvp/pSdhLvPFhbm/msuaA8c5R1F+XmknncpSe161Xhbo0eP5je/+Q3vv/8+vXrVfHtKBYMmTZpw9913s7111dreuJPcezirp7zA8uXL6dGjhxfShS5HtMBRSikni4itRduLb6Nuyy413lbt2rUZOXIkU6ZMoaDA/R94SoWKnj17MmHCBMJcNT92dVa3AYS5Injvvfe8kCy06SBRKaXcOHLkCHvXLKaowP0R/TNx++23c80113D8+HGvbVOpQJWRkcGaNWu8tr3I2Fo0H3A1nTt39to2Q5WeblZKKTc+/fRTlk58lL4PvEqd5p2qvZ3TpxKLg9ZXc8O/Vp+2jk4jpkLRuHHj2LZtG5s3b/baNttfMpYxY7SeakqPJCqllBtTpkwhpm4jaqd29Op2TVEhB35cRWF+nle3q1Qg2bt3L/PmzePaa689Neeytxw+fJjvv//eq9sMNTpIVEqpSmRlZTFnzhwadx/k9V9g+zdm8M3f7yFr0zKvblepQPLxxx9TVFTE6NGjvb7tu+++mxEjRlBYWPmNo8o9HSQqpVQlZsyYQUFBAWd1G+D1bddrdTbhUbHsWbnA69tWKlBMnTqVDh060LGjd4/UA4waNYr9+/ezePFir287VOggUSmlKrFgwQKSk5NJbNrG69t2RUTSoOO57F2zGFOkRzpU6MnOzmbp0qWMGjXKJ9sfPnw4UVFRfPLJJz7ZfigImRtXPPXU9LTcU5/EvILKexGC+16He46ccPveHw8fc7t87Z4ct8uP51beZiMp3n0fxJb1o90ub5EY53Z5o4TK3x/uoUejK8xDn0QPp//cfce8fOZQBam33nqLn3/+mbs//ckn22/UtT+7ls/lwNbV1GvVzSf7UMqpYmNj2b17N7m5uT7Zfq1atRgyZAjTpk3jpZde8volI6FAjyQqpVQlwsLCaNq0qc+236BDb8LCI9m7Rk+HqdAUHx9PvXr1fLb9yy+/nB07drB69WrPK6tydJColFIVePzxx3n44Yd9uo/wqFj6PvBP2o24w6f7Ucpp8vLyGDhwIDNnzvTpfi6//HLWrVtHly41b4QfinSQqJRFIjJVRC4WEa1FBykqKuLNN99k27ZtPt9XYtO2hIW7n95SVZ/WmDMtXLiQ+fPn+/zO48TERNq3b+/TfQQzLRql7HoVuA7YLCLPi0g724EULFu2jN27d3PJJZf4fF9FhQWsn/4aP2fM8fm+QpTWmANNnz6d6OhoBg8e7PN9rVmzhhtvvJG9e/f6fF/BRgeJSllkjEk3xlwPdAe2A3NE5BsRuUVE9PCSJZ9//jlhYWFcdNFFPt9XmCucvWsWs/O7WT7fVyjSGnMeYwwzZsxg8ODBxMbG+nx/eXl5vPfee8yePdvn+wo2IXN3s1JOJSL1gBuAG4EVwPvAecDNwAX2koWuOXPm0LNnT59eUF9agw7nsG3Bxxw7doz4+Hi/7DOUaI3ZV3payuP7M9m+fTvxPS8vM12lb3Tr1o1GjRoxa9YsbrrpJp/vL5jokUSlLBKRacDXQCww0hhziTHmQ2PMvUClowURGSYiG0Vki4g8WsHyUSKySkR+EJEMEdFJTKvIGEOXLl18MgNEZRp0OBdTWMDcuXP9ts9QoTXmPIV5J2jQ4Rzqt03zy/7CwsIYPnw4s2fPpqCg8pZwqryQOZJo3HbN86yg0P37TxS4v/j2cE7lfRI99UH8+sfDbpdv/Omg2+W//JJd6bKEhCi3722bWtft8oFt3X9dYsMr/y8WFe5y+94ID30Uw13u3x8gHbEmGWNOO88oIlHGmFxjTIU/QUXEBUwAhgCZwFIRmW6MWVdqtbnAdGOMEZEuwEeAXotVBSLCa6+95td91m3RmfCoWGbNmuWzxsIhTGvMYRKatKLX2Bf8us+LLrqIt99+myVLlnDeeTqer6qQGSQq5VDPAmUvRvuW4uunKtML2GKM2QogIlOAUcCpX2DGmNJ/ecThvre4KmXv3r00aNDAr413w8IjOKvbBURGum9ur6pFa8xBTFEhecePEFWrjs/3VfpUdn5OPPENm/Hgu4tpuOJ/68y4VweM7uggUSkLRKQR0ASIEZGz+d+BzwSKT4u50wTYWep5JtC7gn1cBjwHNAAurmnmUDFgwAA6d+7Mhx9+6Nf9dr1+HOP1F5bXaI050+HMLSx68TZ63v48DTv39dt+I2LiueD37/ltf8FCB4lK2TEU+BWQDLxU6vWjwGMe3lvRIa5yRzGMMZ8An4hIP+AZoMJeEyJyB3AHQEpKiqfcQaf00YacQ/tZv349tBnglwvqK3LixAmio91Ph6mqxBE1Fur1VVbWxgwAajez07vQFBVijCHMpcOfqtCvklIWGGMmA5NF5ApjzNQzfHsmUHquuGRgl5t9LRSRliJS3xiTVcHyicBEgLS0tJA+ZfbLpmUAfrugvqyLL74YEeHzzz+3sv9g4pQa0/o6XdamDGo1bkFUgvvr3X3h2L6dLH7pTjpf/RCNuw/0+/4DkQ4SlbJARG4wxrwHpIrIg2WXG2NequBtJy0FWotIc+BnYDTFzYJLb78V8GPJRfXdgUjgF699AkEqa9MyIuMSSWjc0sr+mzVrxr/+9S/y8/OJiNAWfjWhNeY8hXm5HPhxFc3Ou9TK/mPrnYUxhqxNGTpIrCJtgaOUHXEl/8YDtSp4VMoYUwD8GpgNrAc+MsasFZGxIjK2ZLUrgDUi8gPFd2leY4wJ+aMY7hhj2L8xg/pteyBhdn40DhkyhOPHj7NkyRIr+w8yWmMOc3DbaooK8qwdqQ9zhVOv9dns35CBfquqRo8kKmWBMeb1kn+fqub7Z1Hmjk1jzGulPv4z8OeaZAw5xtDl6oeIiEuwFmHAgAGEhYWRnp7O+eefby1HMNAac55aZ7Wg8zW/o16rrtYyJLVJY++qr8nO2kVcUhNrOQKFDhJLFHn4q8LT8rz8IrfLD52ovE/ipv0n3L53zZZyl5Gd/v41O9wuLzxaeR/FfYnuZ5Tw9MdW03oxbpc3S8irdFntQven04o87NvT98RdBxOx3EVRRP7hbrkx5j5/ZVHFJCzMr3dbVqR27dr07NmTOXPm8NRT1RrbqBJaY84TlVCXZn19Px+6O/XbFR/FzNq4VAeJVaCDRKXsWGY7gDrdruVziWuQQmJya6s5HnroIU6ccP+Ho6oSrTEHyc8+yu6VC2nUuS+R8bWt5YhLakqb4WOo3ayDtQyBRAeJSllQcuelcoiiwgJW/fsvNE4bQpdrHrKa5corr7S6/2ChNeYsWZuXs+rfzxPfcAJ1LQ4SRYQ2w2+xtv9Ao4NEpSwQkZeNMfeLyAwq7r9m95xMiDm8YwMFudnUb9PDdhQANm/ezJ49e/S6xBrQGnOWrI0ZuKJiHHEEr6iwgINbVxNT7yzbURxPB4lK2fFuyb8vWk2hANi/MQNEqN/G3Uxt/nPfffexffv24sbeqrq0xhxk/8YM6rU62xFNrPOzj/Lt+PtoO+J2QI/cu6MtcJSywBizrOTfBRTPI3sQOAB8W/Ka8qOsjRkkJrchMi7RdhSguBXOhg0byMzMtB0lYGmNOcf27dvJ3p9prfVNWVG16pDQpBVZG/WyVU90kKiURSJyMfAj8A/gFWCLiAy3myq0FBXkc+TnLY451QwweHDx7G5z5syxnCTwaY3Z9/333wOQ5JBBIhTPqnRw62qys7NtR3E0+8d9naKGfTULPbRjKTCVt8g5cqLA7XsPHsxxv+9jh90uZ8+PlS7KD3O5fevx4w3cLj94vPLWPgB5hZV/3oUeetx4bHYaHL1Q/woMMMZsARCRlsBM4AurqUJIWHgEQ579jML8XNtRTuncuTMNGjRg7ty53HKLXmRfQ1pjll199dW8sSmSqAT3Ldf8qX7bNLbOm8KiRYu48MILbcdxLKuDRBEZBvwdcAGTjDHP28yjlAX7Tv7yKrEV2GcrTKhyRUbhioyyHYOR4xed+jg8uTMffz6bQ//4Gilp+jnj3vNsRQtkWmMOEJ1Y33aE09Rt0QVxhTNv3jwdJLphbZAoIi6KpzIaQvFk6ktFZLoxZp2tTEr5i4hcXvLhWhGZBXxE8bHRqyieN1b5yYp/PUP9tj1o2vsi21FO027knYRHx54aIKozozXmDGvXruWxxx7jWOeriW/YzHacU8KjYjj/oTd45pnrbUdxtCoNEkUkGrgbOI/iIlsEvGqMqUnH117AFmPM1pJ9TAFGATpIVKFgZKmP9wL9Sz7eD9Txf5zQtHfvXn7O+JJajVJtRyknVttz1JTWmAPMnj2b6dOnM6jHzbajlJPQpBUREe5n/gp1VT2S+C/gKDC+5Pm1FLcXuKoG+24C7Cz1PBPoXXYlEbkDuAOgaUpKDXanlHMYY/RCMweYN28egGPuuixr55KZ5Bzcp81/q0FrzBnS09Np164dMXXcX99uQ37OMX77299y0UUXMWjQINtxHKmqg8S2xpjSM3LPF5GVNdx3RedQKmp4OhGYCNCjR1pw3KqgVImSo/S3Ah2B6JOvG2PGWAsVQtLT04mIiSexaRvbUSp0cPs6di1Lp9WFNzqiv1wg0hqzJy8vjwULFjBmzBi22w5TAVdkNG+88QbZ2dk6SKxEVVvgrBCRc04+EZHewOIa7jsTaFrqeTKwq4bbVCrQvAs0AoYCCyiug6NWE4UIYwxz5syhXpseiIe7/G2p3zaNgtxsDu/YYDtKINMas2TJkiVkZ2efaunkNGGucPr37096errtKI5V1UFib+AbEdkuItspbkzaX0RWi8iqau57KdBaRJqLSCQwGphezW0pFahaGWMeB46XzDV7MdDZcqaQcPz4cbp160bDTn1sR6lU/dZnA5C1SZv+1oDWmCW5ubn06tWL/v37e17ZksGDB7NlyxZ++ukn21EcqarnL4Z5e8fGmAIR+TUwm+IWOG8ZY9Z6ez8nSYVnt0uvULMz2eFh7sfbYW7uUEyIdv9tqFcv1u3yI0nur/XILsirdFls/SS3742Li3S7vHas+4t+o1yVH6Fx9zWBqnzP3C8OECcbTR4SkU7AHiDVXpzQER8fz/Tp009rO+M0kfG1SUhuTdbGZbQe6rwL/wOE1pglQ4YMYciQIbZjuHXyNPPcuXMZM0avQCirSoNEY4xPhtjGmFnALF9sW6kAMVFE6gCPU3wkPb7kY+VjR44cISEhwXYMjxp0OIdDP6333FxeVUZrzIK8vDyKioqIjo72vLJFHTt2pFOnThw7dsx2FEfSK6GVssgYM6nkwwVAC5tZQklhYSHNmzfnrrvugoZeP1HiVe1G3GE7QkDTGrNj9uzZXH311Xz77bd069bNdpxKiQirV6+2HcOxdO5mpSwSkXoiMl5ElovIMhF5WUScM3dVkFqxYgUHDhygY8eOtqNUmSkqtB0hIGmN2ZGeno6I0L59e9tRqqywUGusLB0kKmXXFIqnCLsCuBLIAj60migEnLybceDAgZaTVM2aqX9n8d/uth0jUGmNWZCenk6/fv2IirI/3aUnR44coXXr1owfP97zyiFGB4lK2VXXGPOMMWZbyeNZoLbtUMEuPT2dLl260LBhQ9tRqiQyLpFDO9aTlZVlO0og0hrzs127drFu3TrHtr4pKyEhARFh7ty5tqM4jg4SlbJrvoiMFpGwksfVwEzboYJZTk4OixYtCqjmufXb9ABjmD9/vu0ogUhrzM9ODrYCZZAIxXc5f/XVV+Tn53teOYToIFEpC0TkqIgcAe4EPgDySh5TgAdsZgt2xhgmTJjADTfcYDtKldVu1p7wqFg90nEGtMbs6d27N88//zxdunSxHaXKBg8ezLFjx1i6dKntKI6idzeX8NSTzxXmfnlEuPvlSTGVX5fRsp77azYyU9zPRZ+X5/5i28Lkyq/Rjolx3+ewtYd9t0ly394gMary7UdFeOgt6eFPGE/fM499Fi0yxtSynSFUxcbGcuutt9qOcUbCXOHUbdVNZ4Y4A1pj9rRp04ZHHnnEdowzMmDAAESE9PR0+vRxboN9f9NBolKWicglQL+Sp18ZYz6vwnuGAX+nuBH9JGPM82WWXw+c/Cl9DLjLGFPT+daDwtSpU+nVqxdNmzb1vLKDpJw7gmvaRlBYWIjLTZN6VZ7WmG+Vbkifc2g/R3Zuon7bHrgind0jsbS6devy+OOPc+6559qO4ih6ulkpi0TkeeA3wLqSx29KXnP3HhcwARgOdACuFZEOZVbbBvQ3xnQBngEmejt7IPrll1+46qqrePvtt21HOWONupzPQw89pAPEM6Q15l97Vi5g6RuPknv0gO0oZ+ypp55y/Awx/qZHEpWy6yKgmzGmCEBEJgMrgEfdvKcXsMUYs7XkPVOAURT/AgTAGPNNqfWXAMlezh2Q5s6dizGGCy+80HaUavnll1/YsmULvXv3th0lkGiN+dH+Dd8Tm5RMbL3GtqOcMWMM69evJzY2ltTUVNtxHEEHiUrZVxs4+Wd3YhXWbwLsLPU8E3A3argV+KJayYLMnDlzSExMJC0tzXaUarnvvvuYO3cuu3fvRjzMfa5OUxutMZ8rKsjnl80/kNzL2bMYlVb6VHlhXi6zHxlOav8r6HDpPaetN+Pe8/wdzRF0kKiUXX8CVojIfEAovm5qnIf3VDQ6qHBiXxEZQPEvsEp/wonIHcAdACkpKVWIHFhO/hIwxjBv6gwSU7tx2atLLKeqnsGDB/PBBx+wdu1aOnXqZDtOoLBaY8FeX6Ud3LaGwrwcktoF5h9hrsgo6rToTNbGZbajOIZek6iUJSISBhQB5wDTSh7nGmOmeHhrJlD6rotkYFcF2+8CTAJGGWN+qWxjxpiJxpg0Y0xaUlLSGX4WgSPnwB5yDu4lqW1g/gIDTvV21Lucq8YJNRYq9QXwy48rkTAX9Vp3tx2l2uq36cGRnzeTe/Sg7SiOoINEpSwpuUbq18aY3caY6caYz4wxe6rw1qVAaxFpLiKRwGhgeukVRCSF4l+INxpjNnk9fACKrXcWQ579jMY9AqfBb1kpKSm0bt1aB4lVpDXmX60vvIn+j71LREy87SjVVr9tDwB+2bzCchJn0NPNJTy0QSTMw/U/UeHux9sJsZX3C+xQL8Htez31aEyuG+N2+aHsyjvIJ3jok9jWQx/ElonufxjER1f+XyzaQ5/EcA+NEj19zwLEHBF5iOK5ZI+ffNEYU+mtgcaYAhH5NTCb4vYcbxlj1orI2JLlrwF/AOoB/yy5dq3AGBO4h9C8JCqhru0INTZo0CDee+898vPziYhwX78K0BrzGwkLI75BYLWWKiuxaVvCo+PI2pRB4+6BMbe7L+kgUSm7xlB8rdPdZV5v4e5NxphZwKwyr71W6uPbgNu8lDHgFRUWsPydJ0k9/7LiKe4C2EMPPcQDDzxAeLj++K4irTE/2Lf+O3avmE/7UXcTGef+wIeThbnC6X3Xi8Q1bGY7iiPo6Wal7OpAcT+2lcAPwHigo81AwejQ9nXsWbmA/JxjtqPUWMuWLWnTpo3e3Vx1WmN+sGflAnb/8BXh0bG2o9RYneadiIzVCXtAB4lK2TYZaA/8g+JfXu1LXlNetH/D9yBh1A/gC+pLmzVrFs8++6ztGIFCa8zHjDHs37CUem26E+YK/CPcRQX5bJnzHvvWBWYXBG/SQaJSdrU1xtxmjJlf8rgDaGs7VLDZvzGD2s3aExEkRwcWLFjA008/zfHjxz2vrLTGfOz4/kxyDuwhqW1P21G8QlzhbFvwMZnfh3zrSx0kKmXZChE55+QTEekNLLaYJ+jkZR/l0E/rSWoXHL/AAIYMGUJ+fj7z58+3HSUQaI35WNbGDICgqTERIaldT/ZvyMAUFdqOY5UOEpWyqzfwjYhsF5HtwLdAfxFZLSKr7EYLDnlHD1A7pR1J7XrZjuI1559/PnFxccyaNcvzykprzMdEhLotuxJbv4ntKF7ToMM55Gcf4dBP621HsSrwLx6oIs/XeLtfIcJVs/F0UYW9+k9y32YmNsL9t6llYuUtbgCO5xdUusxTe5260ZFulye6ae0DUMttCxyX2/d6yubhW1aF77kjBM78VQEqvmEzzvvt67ZjeFVUVBSDBw9m5syZGGP0Jhb3tMZ8rNl5l9LsvEttx/Cq+u16IWEu9q1bQp3moTu7UcgMEpVyImPMT7YzBLOioiIKcrMJjwr8Oy7Luvjii9m4cSNZWVkE+0weNaE15lvHjh3DFBUhHvraBprI2FrUa9WNvONHbEexSgeJSqmgtWzZMr58dAQ9x/4loKfjO+nkPNQApqgtre5+gzFTNgIbT70+495Kp+lWyusee+wx5n8wlQGPf4CEuT87FGh63/1S0A1+z1Rof/ZKqaA2c+ZMiooKSGzSynYUrzv5y8sUFVlOokKVMYaZM2cS36hZ0A0QQWsMdJColApiM2fOpE6zjkTG17YdxSd+XpbOnN9fQl72UdtRVAjauHEjW7dupUGHc21H8Zllbz/B8neesB3DGh0kKqWC0p49e8jIyKBBx+D9BRZTpwF5xw+TtWGp7SgqBM2cOROAhh37WE7iOxGxtdi3/jvy8vJsR7FCB4lKqaD0xRfFjXCDeZBYJ7UjEbEJ7Fv3re0oKgTNnDmTTp06EVO3oe0oPtOww7kU5ubw9ddf245ihQ4SlVJBqV+/fvz1r38lIQivRzxJwlwkte/FvnVLQvq6KWXHI488EvTTQ9Zr052w8MiQ7UmqdzdXkacbnMI9NO2LcdMT0OWhx1lUuPudJ8S4/zYWmsqbNHrad2SE+317yuauF2K4y/2+PbdJ1N5wqnItW7bkwQcfZH6pO4KDUcOOfdi1LJ1DP60L6X5uyv+GDh0KwKQgrrHwqBjqtT6b6dOn8+KLL4ZcT1I9kqiUCjorVqzg008/JT/ffaP5YNCgwzm0GDiayPg6tqOoEDJlyhRWrQqNCWtSz7+Me++9l8LC0JuiTweJSqmgM2HCBG6++WaMm6PowSIithYdLr2HuKTgmRJNOduJEye4/fbbmTBhgu0oftGwU1/uu+8+wsND7+SrlUGiiLwgIhtEZJWIfCIitW3kUEoFn4KCAj777DNGjBhBZKT7aSWDRVFhAVmblpH9y27bUVQISE9P59ixY1x22WW2o/jNgQMHmDp1qu0YfmfrSOIcoJMxpguwCRhnKYdSKsgsWrSIrKwsLr/8cttR/CY/+yhLJjzIziUzbUdRIWDatGkkJCQwcOBA21H85t133+XKK69k8+bNtqP4lZVBojHmS2NMQcnTJUCyjRxKqeAzbdo0oqOjGTZsmO0ofhNVqw71WnZlz6qFtqOoIFdQUMD06dMZOXJkyBypB04dNZ02bZrlJP7lhGsSxwBfVLZQRO4QkQwRydiftd+PsZRSgWjFihUMHTqUuLg421H8qlHXfhzdvY2NGzd6XlmpatqwYQM5OTkhdaoZICUlhZ49e+og0VtEJF1E1lTwGFVqnd8DBcD7lW3HGDPRGJNmjElLqp/kq7hKqSCxcOFCJk+ebDuG3zXq0g+ATz75xHISFcw6depEVlYWI0eOtB3F7y6//HK+//57du7caTuK3/jsVh1jzGB3y0XkZmAEMMg44BZEj62PjPsVPLT8c/uVDvPQhNFTP8EiD18+d4vDPHzinr4uES732V1umh166j3pqQ9iiLWrUlVgjEFESExMtB3F72LqNKB2sw7MnTuXRx991HYcFYRO1ldMTIztKFZcfvnljBs3jvnz53PTTTfZjuMXtu5uHgY8AlxijMm2kUEpFVwKCgro2LEjb7zxhu0o1vQY8/Sp6QiV8rYvvviCrl278uOPP9qOYkWbNm3YunVryAwQwd6MK68AUcCcku7lS4wxYy1lUUoFqJGlZnrYt/471q9fz+vf7Wf6ieCdAcKdmDoNQ7KXm/KNkWVmUlk++WX2b9nOfdN3EBYemu2WmjdvbjuCX9m6u7mVMaapMaZbyUMHiEqdAREZJiIbRWSLiJQ7tygi7UTkWxHJFZGHbGT0t13L0gmPiadBh3NsR7Fq8uTJnH/++RTpXM41ojV2uoLcHPau/pqzzh5AWHiE7TjWFBUVcd111wX9nNUnOeHuZqXUGRARFzABGA50AK4VkQ5lVjsA3Ae86Od4VhTmnWD3ygWc1bU/rojQactRkfDwcBYtWsTixYttRwlYWmPl7V2zmMK8EzTp4fZ2g6AXFhbGgQMHePPNN0PiDzE9L6FU4OkFbDHGbAUQkSnAKGDdyRWMMfuAfSJysZ2I/rV3zWIKc3Noknah7SjWXXrppcTFxfHee+9x/vnn244TqLTGyvg540ui6zSgbosutqNYUfrU+94GPdk+ezbnP/gadVv+7+sx497zbETzKT2SqFTgaQKU7sGQWfJatZzWi3R/YPYijW/UnBYDrqFeq662o1gXFxfHZZddxkcffcSJEydsxwlUXquxYKgvgCZpF9Jm2BjEU1uKENCoy/m4IqPJXDrbdhSf0++2UoGnouY/1W4jdVov0qTA7EWa0LgFHS77NRLmsh3FEX71q19x6NAhPv74Y9tRApXXaiwY6gugSY/BpJwbEgdNPQqPiuWsbhewa1k6BbnB3aBFTzdXkeeefNXvo+ipV6ErzMPPJk8/utxs3lMvQk/ctEEs3r6bz037HFZbJtC01PNkYJelLNbtWfU1UYn1qNOs7CVjoWvAgAHcc889tGnTxnaUQKU1VsIUFbH962k07j6IqFp1bMdxjNTzLye6dhJFhYW2o/iUDhKVCjxLgdYi0hz4GRgNXGc3kh1FhQWs/s9LJDZpRa+xL9iO4xhhYWG88sortmMEMq2xElkbM1g79e9ExtcO+ZtWSqvdrD21m7W3HcPn9HSzUgHGGFMA/BqYDawHPjLGrBWRsSIyFkBEGolIJvAg8H8ikikiCfZS+8b+dd+ReziLlD6hN0VYVWzcuJHPPvvMdoyAozX2Pzu+nUFEbMKpaR/V/5iiQvauWcyxvT/ZjuIzeiRRqQBkjJkFzCrz2mulPt5D8SmyoLZtwX+ITkyiQcc+tqM4Rum7MJe/8wT713/PoKenER71v6nUgvEuTG/TGoOcA3vZs+prml9wVci3lqpIwYlslr39BE16DKbrdcE5FaYeSVRKBaRVq1aRtWkZqf2uIMylf+9WJLXfFeTnHGPnd7M8r6xUGdsWFt/41LzflZaTOFNEbC2Sew3j54w55B45YDuOT+ggUSkVkDZs2EBUYn1S+l5iO4pj1WnemdqpHdk2/yNMUXBfYK+878ThXzir2wXE1G1oO4pjtbjgaooK8ti+6BPbUXxCB4lKqYB09dVXM+jJ/xAZW8t2FMcSEVoMuIbsX3axZ3Vozmetqq/7zX+g243/ZzuGo8U3TKFhp7789PUnZGcHXzscHSQqpQLOhg0bMMboaeYqOKtrP+IbpZKdFZIdXFQ15OXlsWXLFgCtsSpoMXA0rqgYfvzxR9tRvE4HiV4i4v4RJlLthyvMw8Pl4eHmvWFhuH142reIp0flXxOlqiMrK4u0tDQee+wx21ECgoS56PfI27QcdK3tKCpATJ48mXbt2nHk5y22owSEui27MuDxf9O5c2fbUbxOB4lKqYDy4osvkp2dzU033WQ7SsA4eTTo0I7iI7BKVSYvL49nn32WtLQ0ajVuaTtOQBARwlzhnDhxglWrVtmO41U6SFRKBYz9+/fzyiuvcO2119K+ffA3svWmvasXs+jF29m37lvbUZSDvfPOO+zYsYMnn3zS7YxZqrybbrqJYcOGkZOTYzuK1+ggUSkVMJ5++mlycnJ4/PHHbUcJOEkdehNbrzEbZ06iqKjIdhzlQMeOHePJJ5/knHPOYejQobbjBJy7776b3bt3889//tN2FK/RQaJSKiCcOHGC2bNnM3bsWNq1a2c7TsAJc4XTdsTtHMnczFtvvWU7jnKgb775hgMHDvDSSy/pUcRquOCCCxg2bBhPP/00e/futR3HK3SQqJQKCNHR0axatYrnnnvOdpSA1bj7IOq27Mq4ceM4ePCg7TjKYS688EJ27tzJueeeaztKwHr55ZfJyclh3LhxtqN4hd7brpRyrJNTzB3euYn4him4IqMtJwpsIkLHK37DhrcfYe3atZx3nk7Pp8AYwzfffEPfvn1JSkqyHSegtW3blvvvv5/FixeTm5tLVFSU7Ug1ooNEB/B0VF9wv4KnmxX1rIEKZLlHD/Ldq7+lXqtu9BjzjO04AS8xuTU7duwgOloH3KGs9BzfmUu/5Id3n6Hn7c/TsHNfi6mCwzPPPENERARhYYF/slYHiUopR1s79WXyc47RetivbEcJGtHR0RQVFTFp0iSuvfZaatXSWWtC1Ykjv7B26svUSe1Eg47n2I4T0EoPvKH4D9z9678judewU6/NuDewjt4H/jBXKRW0Mr//L7uWz6PNsFtI0J5tXrVy5UrGjh3LAw88YDuKssQUFbLyvT9RmJ9L1+vHIWEu25GCyo9zP+CH9/9E1qbltqNUmw4SlVKOtHLlSlZ9+CJ1W3Wj5eDrbMcJOmeffTaPPvoob775JpMmTbIdR1mw6Yu32b/hezpefh/xDVNsxwk6bYbfQlxSMsvfeYKcg4F5t7OeblZKOVLt2rVJatuTLqN/p/PHetnJ02Km4VCS2qVzx9i7eGddAXWadwIC75SYqp5ajVuS2u8KUvpcYjtKUAqPiiXttj+x+K93kvHGY/S5P/D6J+qRRKWUoxw/fpyCggKaNWtGzzueIyqhru1IQUvCXJz9q6eIqdOQZW/9gcL8XNuRlB8cPnwYgMZnD6DTlfdrT0QfqtUolbNv+gOHf97Mhs9ftx3njOkgUSnlGNnZ2VxyySVce+21Osewn0TG1qL3XS/S9fpxuCICu12H8iwjI4MWLVqwe+UC21FCRsPOfel+85O0Hnqz7ShnTAeJSilHyMrKYtCgQcyfP59Ro0bp0Q0/iktKJqldT6D4ZqHlywP3QntVuS+//JIBAwZQq1Ytaqfo3Of+1Lj7QCLjEsnLy+Opp57i+PHjtiNViQ4Sg4CI+4dSTrd27Vr69u3LihUrmDp1KjfccIPtSCGpMD+XTV+8Tb9+/fj0009tx1FeYozh9ddf5+KLL6ZFixZ88803xNRpYDtWSFq4cCFPP/00F1xwAZmZmbbjeKSDRKWUFSPHL2Lk+EVc/PJX9Ow3hO2799Nj7Eu8lZl0apnyL1dEFH3un0C7du247LLLuOeee8jJybEdS9XQV199xdixYxk0aBALFy6kcePGtiOFrMGDB/Ppp5+yfv16unTpwrRp02xHcktvGVRKWXF09zbikpIJC4+g+81PEFOvEdEJ9WzHCnnRifVZvHgxjz32GC+99BLz5s1j5cqVREZG2o6mqmjk+EWYoiKO7t5KQpNWGBNO2u3PEd6xDzf8a7XteCFv5MiRLF++nOuvv54rrriCJ554gieffNJ2rArpIFEp5Vc//PADzz33HAv+8x86jLqbFgNHU6d5R9uxVClRUVH89a9/ZcSIEacNEOfPn0///v2DYrqxYFVQUEDm0i/5Mf09jmf9zIDff0BM3YY06qxtjZyg9BmSBte/QNtGHzAvuynLxi8i5+BeCvNziW+Q4pg2VFYrXUQeEhEjIvVt5lAq0IjIMBHZKCJbROTRCpaLiPyjZPkqEeluI+dJxhhee+01+vTpw9lnn80XX3xBy0HXkdx7uM1YqhInT/e/tCaCua40Ro5fxHkPvsbAgQOp1TCFNsNvYcOGDUF9B3qg1diuXbv43e9+R7Nmzfjh3eI5zrteN47o2kk2Yyk3wlzhtL7wJhKbtgVgy5fv8tUfb2DJhAd55513OHTokN2AWBwkikhTYAiww1YGpQKRiLiACcBwoANwrYh0KLPacKB1yeMO4FV/5SssLGTjxo18+OGHjB8//mRm/v3vf3P06FH+8pe/sGPHDtpfMpbIuER/xVI1lJjSjrNvfoKYug3ZPHsy7du3p3nz5qxfvx6AQ4cOceLECcspvcPpNXb8+HG+/fZbXn31VebPnw9AXl4e//jHP0hLS6Pn7c/T75F3aNJjMKJHfQNGm+FjaD30ZrKzMrnllluoX78+l1566anl+/bt8/sfZjZPN/8NeBj4zGIGpQJRL2CLMWYrgIhMAUYB60qtMwr4lyn+ibJERGqLyFnGmN3uNmyMITs7m/z8/FOPpKQkwsPDycrKIjMzk8OHD3Po0KFT/44dO5bIyEj+9re/MXHiRH766adTNzu4ImP4Ir8DrogoYi5+jGbRsSwUYeG7a3zyhVG+E+YKp0mPwTTpMZgTh7O4tM7PfPXVV6SmpgLw5z//mRdeeIHU1FRSUlJo2rQpKSkpPP3004gIK1eu5ODBgyQkJBAVFUVkZCQxMTEkJycDxYMcAJfL5YTT2T6rsby8vNPqS0SoX7/4ZNrmzZs5dOjQafXVqFEjRowYAUBSu54c27eDE4f2Q8lgIbn3cLpdHwHAgGc+oygmnobe/Eoov4lKqEvbi26lzfAxHNq+jn3rvmV9fvipU9Tpj19OwYnjxNQ7i5g6DbmiX1cuvPDCUwPJ//73vyQkJBAbG0tUVBRRUVHUrVuX2rVrY4whNzcXl8uFy1X1ObqtDBJF5BLgZ2PMSu2FptQZawLsLPU8E+hdhXWaAG5/gS1fvpy4uLjTXtu6dSvNmzdn0qRJjBs3rtx7Pj/WnKhaddiZsY9DUQ1p1LsLCY1bkJDchlqNUgkLL/4FFhETV+69KjBFJ9bnv0X1oV9Xrp60DIADOU1pMfgGju/P5Ifte/n2h7UUnMjmhwZDAVg++Wl2LZtz2nYaNGjA3r3Fc9peddVVTJ8+HYArr7zSj59NhXxSY+vXrycq6vSG5X379mXRouLT+1/98QaO7f3ptOVJ7XrSe1ttAKJq1SE6sT6x9RqT0KQVCcmtianzvyFhREx8FT895WQiQp3mHU+7VtsYQ+vhv+Lorq3kHNhDzsF9fPjhhwBceumlFBQUMHx4+ct3fvvb3/Liiy9y9OhREhPP/MyN+OrQpYikA40qWPR74DHgQmPMYRHZDqQZY7Iq2c4dFB/KB+gEOP0QRH2gws/FIZyeDwIjY1tjTC0bOxaRq4ChxpjbSp7fCPQyxtxbap2ZwHPGmEUlz+cCDxtjllWwvUCqsUD4v6EZvSMoaizA6gsC4/+GZqy5KtWXz44kGmMGV/S6iHQGmgMnjyImA8tFpJcxZk8F25kITCx5b4YxJs1Xmb3B6Rmdng8CJ6PF3WcCTUs9TwZ2VWMdILBqzOn5QDN6S7DUWCDVF2hGb3F6xqrWl98v/DDGrDbGNDDGpBpjUikutO4VDRCVUhVaCrQWkeYiEgmMBqaXWWc6cFPJHZjnAIc9XSullDpFa0wptE+iUgHHGFMgIr8GZgMu4C1jzFoRGVuy/DVgFnARsAXIBm6xlVepQKM1plQx64PEkqOJVTXRVzm8yOkZnZ4PNKNHxphZFP+SKv3aa6U+NsA91di007/2Ts8HmtFbgrHG9OvuHZqx5qqUz2c3riillFJKqcBlvRmVUkoppZRynoAdJDp1Sj8ReUFENpRM0/SJiNS2nekkT9NM2SYiTUVkvoisF5G1IvIb25kqIiIuEVkhIp/bzuIrTq0v0BqrrkCpL9Aas82pNebk+oLgrLGAHCSKs6f0mwN0MsZ0ATYB5bsPW1DFaaZsKwB+a4xpD5wD3OPAjAC/AdbbDuErDq8v0BqrrkCpL9Aas81xNRYA9QVBWGMBOUjkf1P6Oe6CSmPMl8aYgpKnSyjuneUEp6aZMsbkASenmXIMY8xuY8zyko+PUvwfuIndVKcTkWTgYmCS7Sw+5Nj6Aq2x6gqE+gKtMSdwaI05ur4gOGss4AaJUmpKP9tZqmAM8IXtECUqm0LKkUQkFTgb+M5ylLJepviHe5HlHD4RYPUFWmPV4uD6Aq0xp3FKjQVMfUHw1Jj1FjgVkSpM6effRKdzl88Y81nJOr+n+NDz+/7M5kZFk2Q78q9YEYkHpgL3G2OO2M5zkoiMAPYZY5aJyAWW41Sb0+sLtMZ8yan1BVpj/hSANRYQ9QXBVWOOHCR6a0o/f+c7SURuBkYAg4xzegxVeZo2m0QkguLiet8YM812njL6ApeIyEVANJAgIu8ZY26wnOuMOL2+QGvMVxxeX6A15jcBWGOOry8IvhoL6D6JIrIdSDPGOGYSbREZBrwE9DfG7Led5yQRCaf4AuRBwM8UTzt1nTFmrdVgpUjxT83JwAFjzP2W47hV8hfYQ8aYEZaj+IwT6wu0xqorkOoLtMZscmKNOb2+IDhrLOCuSQwArwC1gDki8oOIvObpDf5QchHyyWmm1gMfOam4SvQFbgQGlnztfij5a0ep0rTGqkfrS1WV42osAOoLgrDGAvpIolJKKaWU8g09kqiUUkoppcrRQaJSSimllCpHB4lKKaWUUqocHSQqpZRSSqlydJColFJKKaXK0UGiUkoppZQqRweJSimllFKqHB0khggR6Skiq0QkWkTiRGStiHSynUupYKD1pZTviUiqiGwQkckl9faxiMTazhXMtJl2CBGRZymeqzEGyDTGPGc5klJBQ+tLKd8SkVRgG3CeMWaxiLwFrDPGvGg3WfDSQWIIEZFIiue7PAH0McYUWo6kVNDQ+lLKt0oGiQuNMSklzwcC9xljLrWZK5jp6ebQUheIp3hOzmjLWZQKNlpfSvle2SNbeqTLh3SQGFomAo8D7wN/tpxFqWCj9aWU76WIyLklH18LLLIZJtiF2w6g/ENEbgIKjDEfiIgL+EZEBhpj5tnOplSg0/pSym/WAzeLyOvAZuBVy3mCml6TqJRSSinHK7km8XNjjHYO8BM93ayUUkoppcrRI4lKKaWUUqocPZKolFJKKaXK0UGiUkoppZQqRweJSimllFKqHB0kKqWUUkqpcnSQqJRSSimlytFBolJKKaWUKuf/AQFp89dedZeCAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 792x216 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 3)\n",
    "axs[0].hist2d(x_list, v_list, bins=21, range=[[-4, 4], [-4, 4]], cmap=\"Blues\")\n",
    "axs[0].set_xlim(-4, 4)\n",
    "axs[0].set_ylim(-4, 4)\n",
    "axs[0].set_xlabel(\"x\")\n",
    "axs[0].set_ylabel(\"p\")\n",
    "axs[0].set_xticks([-4, -2, 0, 2, 4])\n",
    "axs[0].set_yticks([-4, -2, 0, 2, 4])\n",
    "axs[0].set_aspect(1)\n",
    "\n",
    "position = np.linspace(-4, 4, 500)\n",
    "boltz_factor = np.exp(-0.5/kT*k_spring*position**2)\n",
    "Z = np.trapz(boltz_factor, position)\n",
    "boltz_factor /= Z\n",
    "axs[1].plot(position, boltz_factor, color=\"k\", linestyle=\"--\")\n",
    "axs[2].plot(position, boltz_factor, color=\"k\", linestyle=\"--\")\n",
    "\n",
    "axs[1].hist(x_list, density=True, bins=21, alpha=0.8)\n",
    "axs[1].set_xlim(-4, 4)\n",
    "axs[1].set_ylim(0, 0.5)\n",
    "axs[1].set_box_aspect(1)\n",
    "axs[1].set_xlabel(\"x\")\n",
    "axs[1].set_ylabel(\"probability\")\n",
    "axs[1].set_xticks([-4, -2, 0, 2, 4])\n",
    "\n",
    "axs[2].hist(v_list, density=True, bins=21, alpha=0.8)\n",
    "axs[2].set_xlim(-4, 4)\n",
    "axs[2].set_ylim(0, 0.5) \n",
    "axs[2].set_box_aspect(1)\n",
    "axs[2].set_xlabel(\"p\")\n",
    "axs[2].set_ylabel(\"probability\")\n",
    "axs[2].set_xticks([-4, -2, 0, 2, 4])\n",
    "\n",
    "fig.set_size_inches(11, 3)\n",
    "plt.savefig(\"md_ho_nhc.pdf\", bbox_inches=\"tight\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ovito",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
